#!/bin/bash
#
# Script to derive the lattice constant of a hcp crystal.
# It needs the fitted equilibrium volume data, which can
# be obtained by:
#    eos_fit ev.dat 2 3

vol=22.7601644494892

cat > .plot.script << EOF
f(x) = a * x**3 + b * x**2 + c * x + d

fit f(x) 'ev.dat' u 2:4 via a, b, c, d

a0 = f(${vol})

ca = 2*${vol} / (3**0.5/2*a0**3)

print ""
print "The equilibrium lattice constant of the hcp metal is:"
print "   a   = ", a0, " Angstrom"
print "   c/a = ", ca
EOF

gnuplot .plot.script
rm -rf .plot.script fit.log
